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We study a system of few colloids confined in a small spherical cavity by event driven molecular 
dynamics simulations in the canonical ensemble. The colloidal particles interact through a short 
range square-well potential, which takes into account the basic elements of attraction and excluded- 
volume repulsion of the interaction among colloids. We analyze the structural and thermodynamic 
properties of this few-body confined system in the framework of the theory of inhomogeneous fluids. 

Pair correlation functions and density profiles across the cavity are used to determine the structure 
of the system and the spatial characteristics of its inhomogeneities. Pressure on the walls, internal 
energy and surface quantities such as surface tension and adsorption are also analyzed for the 
whole range of densities, temperatures and number of particles considered. We have characterized 
the structure of systems from 2 to 6 confined particles as function of density and temperature, 
identifying the distinctive qualitative behaviors all over the thermodynamic plane T — p in a few- 
particle equivalence to phase diagrams of macroscopic systems. Applying the extended law of 
corresponding states the square well interaction is mapped to the Asakura-Oosawa model for colloid- 
polymer mixtures. We link explicitly the temperature in the confined square-well fluid to the 
equivalent packing fraction of polymers in the Asakura-Oosawa model. Using this approach we 
study the confined system of few colloids in a colloid-polymer mixture. 


I. INTRODUCTION 

During last decades colloid physics has been one of 
the areas of more intensive activity in the field of soft 
matter. Colloidal dispersions have wide and well known 
applications in the chemical, pharmaceutical, medical 
and food industries.[1] Moreover, colloids posses certain 
properties that make them good model systems for basic 
research.[2-5] The size of colloidal particles, in the range 
lnm — 10 /im, make possible its direct experimental ob¬ 
servation giving rise to beautiful and conclusive experi¬ 
ments. It has enabled the direct study of phase transi¬ 
tions like the gas-liquid, the gas-solid and gelation ones, 
which is impossible for simple fluids. [6-9] 

Some colloidal suspensions behave like hard spheres 
(HS) where the unique relevant feature of the interacting 
potential is the repulsion length a (particle diameter), re¬ 
lated with the excluded volume.[10, 11] Free energy and 
pressure of few HS colloids confined in a spherical pore 
were studied recently by molecular dynamics and theo¬ 
retical approaches. [12-14] In several cases, it is necessary 
to complement this hard-core repulsion with an attrac¬ 
tive, short-range interaction, to properly describe the po¬ 
tential between colloidal particles. For this purpose, the 
square-well (SW) potential has been utilized.[15] Systems 
of particles that interact through SW potential have been 
investigated extensively, exploiting the fact that its sim¬ 
plicity enables to study both, in bulk and in confinement, 
not only through Monte Carlo and molecular dynamics 
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simulations, but also with analytic approaches. SW sys¬ 
tem is the simplest interaction model that includes a re¬ 
pulsive core and an attractive well with tunable range, 
giving rise to phase transitions and coexistence regions. 
The thermodynamic properties, phases coexistence,[16- 
22] structure,[23] crystallization, glassy behavior[24] and 
percolation phenomena[25] of short-range SW fluids were 
extensively studied in bulk. Properties of inhomogeneous 
SW systems at free interphases and in confinement, were 
also studied. [25-29] Molecular dynamics studies to obtain 
bulk free energy of short range SW particles were per¬ 
formed recently. [17, 30] Pressure on the wall and struc¬ 
tural properties were studied for two SW particles in 
a spherical cavity. [31] Liquid-vapor coexistence of SW 
fluid confined in cylindrical pores [32] and short range SW 
potential in the context of effective interactions among 
proteins[33] were also studied. 

Adding long, flexible, non-adsorbing polymer chains to 
a colloidal suspension causes changes on the phase behav¬ 
ior of the colloidal system. We shall focus on the rather 
simple Asakura-Oosawa model (AO) for colloid-polymer 
mixtures. [34, 35] In this model based on HS-type inter¬ 
actions, the colloids are taken as hard spheres with di¬ 
ameter c t. The polymers with diameter cr p are excluded 
by a center of mass distance of a+ ° F from the colloids. 
However, polymers are treated as non-interacting parti¬ 
cles that can overlap. This kind of interactions leads to 
an effective attractive two-body potential between col¬ 
loids, due to an unbalanced osmotic pressure arising from 
depletion. [36] Although, for a small enough diameters ra¬ 
tio (7 p /a the two-body effective potential has a short- 
range attractive well, and the three- and many-body po¬ 
tentials are null. [37] In addition, when cr p /a is small the 
main features of the colloid-colloid effective pair poten- 
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tial are similar to the simpler SW potential. Studies of 
the AO model with different values of diameter ratio have 
shown interesting properties in bulk, as glassy states and 
demixing, [38, 39] as well as in confinement.[40-42] 

In this work, we study thoroughly highly confined col¬ 
loids in spherical pores, that show different demeanor 
as compared to their bulk counterparts. The confining 
cavity (thought as a nanopore) breaks the translational 
symmetry of a system and causes the appearance of spa¬ 
tial variations. Even farther from bulk case, our aim is 
to work with low number of particles and a system size 
comparable to its constituents elements. This implies 
that we are far away from thermodynamic limit. Nev¬ 
ertheless we work in the frame of statistical mechanics. 
In addition, low N systems allow for analytic solutions 
at some degree, thus letting us to make a direct com¬ 
parison between theory and simulation results. On this 
line, we study the system of few short-range colloids in a 
hard-wall spherical pore. The pure colloidal suspension is 
studied by event-driven molecular dynamics simulations 
using the SW model in a wide range of number densi¬ 
ties and relevant temperatures. We also map the simu¬ 
lated system to a colloid-polymer mixture in a spherical 
pore. In this case we adopt the AO model and consider 
the case of small ratio a p /<j. We connect the short-range 
SW potential and the AO model by adopting the effective 
short-range colloid-colloid pair potential through an ex¬ 
tended corresponding-states law, following that of Noro 
and Frenkel.[43, 44] 

The paper is organized as follows, in Sec. II we pro¬ 
vide details of the interaction model and the statistical 
mechanics theoretical grounds for both few body sys¬ 
tems: colloidal suspensions and colloid-polymer mixtures 
in pores. In Sec. Ill we present the simulation technique 
and the way in which we obtain canonical ensemble sim¬ 
ulations at constant temperature. Sec. IV is devoted 
to present the density profiles, pair correlation functions, 
thermodynamics properties and phase diagrams of sys¬ 
tems of 2, 3, 4, 5 and 6 colloidal particles in a spherical 
cavity. There, we analyze the results for the simulated 
SW system and also for the equivalent AO system, when¬ 
ever possible. We present a final discussion and conclu¬ 
sions in Section V. 


II. THEORETICAL BACKGROUND 

We study a system of N colloidal particles of diameter 
<7 confined in a spherical cavity of radius i? D , at constant 
temperature T. The particles interact with the cavity 
through a hard wall potential which prevents them from 
escaping to the outside. Thus, the effective radius of the 
cavity is i? e ff = R 0 — cr/2, which represents the maxi¬ 
mum possible distance between the center of the cavity 
and the center of each colloid. The temperature of the 
fluid is determined by the wall temperature T, which is 
fixed. We adopt the effective volume V = 47ri? e ff/3 to 
measure the size of the available space for particles and 


consistently define the mean number density p = N/V. 

Given that we deal with a small number of colloids, 
the ensembles equivalence does not apply. The statistical 
mechanical and thermodynamic properties of the system 
with constant N and T is obtained from its canonical 
partition function (CPF), Qn- One actually works with 
the configuration integral (Cl), given that kinetic degrees 
of freedom integrate trivially. Thus, the CPF reads 

Qn = -^A~ 3N Z n . (1) 

Here A is the thermal de Broglie wavelength and Zn 
is the Cl of the system. For pair interacting particles 
Zpf = J v Y[ <jk> e jk dr N , where the Boltzmann factor for 
the j,k -pair is ej k = exp [— /3 <fi{rj k )\, rj k is the distance 
between both particles, <j) is the pair potential and the 
inverse temperature is defined as /3 = 1 /kT, with k the 
Boltzmann constant. For a system in stationary condi¬ 
tions with fixed N and T the Helmholtz free energy (F) 
reads 

F = U-TS = -p~ 1 lnQ N , (2) 

where U is the system energy and S its entropy. The 
reversible work done at constant temperature, to change 
the cavity radius between states a and b is 

b 

F b -F a = -J P w dV , (3) 

a 

with dV = A dR e g. Here, P w is the pressure on the spher¬ 
ical wall which is an EOS of the system. The derivative 
of Eq. (3) at constant T gives the pressure on the wall 
through 

P w = -A-'dF/dRee , (4) 

which meets the exact relation known as contact 
theorem [45, 46] 

(3P W = p (R eS ) . (5) 

In this ideal gas-like relation, p (R e g) is the value that 
takes the density profile at contact with the wall. This 
extended version of the contact theorem for planar walls 
applies to curved walls of constant curvature (spheres 
and cylinders), for both open and closed systems. A 
complete discussion of the presented statistical mechani¬ 
cal approach for few-body confined system that includes 
other properties, such as the energy, may be found in 
Refs. [31, 471. 

Here we consider the confined colloids as particles that 
interact through the square well potential: 

! oo if 0 < r < er, 

-£ if a < r < (1 + A) a, (6) 

0 if r > (1 + A) a, 
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where e > 0. In this work, we study the short-range 
SW system with A = 0.1. The bulk and interfacial 
properties of short-range SW fluids have been studied 
elsewhere. [20, 25, 44] For reference, we note that the 
bulk SW system with A = 0.1 has a metastable fluid- 
vapor transition with critical temperature T = 0.47 s/k 
and density p = 0.47 <t _ 3 .[22, 44] 

Statistical mechanics of the semi-grandcanonical 
confined AO system 

In the AO model a particular colloid-polymer mixture 
is characterized by the ratio q = — . The diameter of the 
polymer coil is given by a p = 2 R g with R g the radius of 
gyration of the polymer. The effect of temperature on R g 
was studied in Ref. [48]. In the AO model the tempera¬ 
ture does not play any relevant role and thus we consider 
T as fixed to fix cr p . We consider the mixture of N col¬ 
loidal particles and polymers at chemical potential p p , 
at a given temperature. The spherically confined system 
of colloids is such that the polymers, which are much 
smaller than colloidal particles, can freely pass through 
the semipermeable wall that only constrains the colloids 
into the pore. The colloid-polymer AO mixture is an 
inhomogeneous system that can be analyzed in the semi¬ 
grand canonical ensemble. Its partition function is 

A -3N N p 

Sm = _ /vr^AT ~\ Zn ’ n »' (7) 

' N p p ' 

with z p = A~ 3 e^ fJ,p , A p the thermal de Broglie length of 
the polymer and A the same magnitude for the colloidal 
particle. Zn,n p is the Cl of the mixture with the N 
colloids, whose centers are constrained to the pore with 
volume V = 47rA 3 ff /3, and the N p polymers in the larger 
volume V p . In the Appendix A it is shown that Eq. (7) 
transforms to 

■= _ „-(p P v exc N) ry(AO) 

“m -“p jy, e > Vv 

where p p is the mean number density of the pure (ho¬ 
mogeneous) polymer system and Z^ 0 ' 1 is the Cl of N 
confined colloids interacting trough the effective pair- 
potential </>AO- Vexc stands for the excluded volume de¬ 
fined in the Appendix A. Given two colloids at a distance 
r apart, </>Ao( r ) is infinity for r < a and it is zero for 
r > a(l+q). In the attractive well region a < r < a(l+q) 
it is 


( 9 ) 

where x = r/a and rj p = p p (na p /6) is the packing 
fraction of the polymers (note that p p takes any pos¬ 
itive value). The expression in Eq. (9) is minus p p 
times the volume of intersection of two spheres of radius 
a (1 + q) /2 whose centers are at a distance r. 


For the thermodynamic analysis of the AO model we 
use as reference the pure polymer system (an ideal-gas) 
with the grand-free energy given by tt p = U p — TS p — 
p p N p (with N p = PpVp). The semi-grand free energy 
of the mixture is Q m = U m — TS m — p p N p and given 
that the system is athermal its energy is purely kinetic 
U m = | kTN + | kTN p . We define 

FaO = ^m — Qp > 

= U c -T(S m - S$) + (3fcT/2 - Pp) ANp , (10) 

with U c = | kTN , AN p = N p — N p and O = —/31n H. 
Here, Fao is the free energy of the confined colloids in 
the polymer solution as an excess over the pure polymer 
system. The pressure exerted by the colloids on the wall, 
the osmotic pressure, is 

P w = ~ A^ 1 d,F ao /dR e ?i , ( 11 ) 

and relates with the colloids density distribution p (r) 
through the contact theorem 

pP w =p(R e s). (12) 

Eqs. (11, 12) are essentially the same that Eqs. (4, 5) but 
the meaning of each magnitude corresponds to different 
systems. 

Extended law of corresponding states 

The short range SW potential and its capability to de¬ 
scribe any short range potential ( universality ) was pro¬ 
posed by Noro and Frenkel in his extended version of 
the corresponding state law. [43] It is based on a map¬ 
ping between different systems using three parameters: 
the effective hard core diameter, the well depth and the 
adimensional second virial coefficient. The later was pro¬ 
posed as a measure for the range of the attractive part 
of the potential. The scheme was used previously for 
studying the critical properties of the liquid-vapor transi¬ 
tion for interaction models including Lennard-Jones and 
Hard-Yukawa.[43, 49, 50]. Recently, it was employed to 
analyze the behavior of proteins in water solutions.[44] 

It is important to note that we apply the law of corre¬ 
sponding states to analyze confined systems composed 
by few particles. This is very unusual and thus we 
made some checks to validate the overall approach that 
will be shown in Sec. IV. Our application of the law 
for hard-core systems is based on the use of two natu¬ 
ral scales: the hard core diameter for the length scale 
and the depth of the attractive well for the temperature 
scale. The reduced second virial coefficient is given by 
B = — 257 / [exp (— /3<f>(r)) — l] dr with b 2 = 2 tt<j 3 /3. 
For the SW system, this gives explicitly: 

Rsw = 1 - (e 1/T * - l) 3A (1 + A + A 2 /3) , (13) 

where we have introduced the adimensional temperature 
T* = Tk/e (its inverse is /3* = 1/T*). In the limit of a 


3a; a: 3 

2(1 + q) + 2(1 + q) 3 


P<f>Ao{x) =-rj p (l + q 1 ) 3 1 
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q = 2A 

Cl 

C2 

C3 

0.2 

0.886964 

0.22724 

0.015641 

0.15 

0.884047 

0.29956 

0.013871 

0.1 

0.876335 

0.44847 

0.012614 

0.05 

0.847353 

0.92875 

0.013498 


Table I. Fitting parameters for the SW-AO mapping. They 
relate T* of the SW system to the value of the packing frac¬ 
tion r] p (see Section II) in the corresponding state of the AO 
system. 


T 

T* (SW, A = 0.1) 

r, p (AO, q = 0 . 2 ) 

0.0051 

0.2 

0.81071 

0.0675 

0.4 

0.43466 

0.1758 

0.6 

0.29597 

0.3033 

0.8 

0.22369 

0.4395 

1.0 

0.17939 

0.5805 

1.2 

0.14950 


Table II. Law of corresponding states for short range poten¬ 
tials of SW and AO types. First column presents the sticky- 
sphere temperature parameter r. Second and third columns 
show the temperature for the SW particles and the corre¬ 
sponding packing fraction of polymer for the AO model, re¬ 
spectively. 

very narrow well, the SW potential reaches the Baxter’s 
sticky spheres limit with 

B = 1 - 1/4 r. 

Here r (that grows monotonically with T) plays the role 
of temperature. [51, 52] In the AO model the effective 
colloid-colloid second virial coefficient Bao is not ana¬ 
lytically integrable. To relate B sw with Bao, we link the 
width of the wells by q = 2A. This mimics the fact that 
P4>ao (x) is deeper near the hard-core of the particle. For 
the AO system with corresponding well-range q = 2A, the 
relation Bao (VpiQ = 2A) = B sw (T*, A) fixes the pair of 
equivalent states T* <—> rj. This mapping allows the nu¬ 
merical evaluation of Bao and the fit of a nearly linear 
relation between 1/T* and 77 for several values of A. The 
coefficients obtained for the fitted function 

1 2 3 

= 3 {l + l-Sq^^Cihi, (14) 

i= 1 

with hi = r] p , h .2 = rjp and /13 = ?7p/(0.1 + rj p )^ are shown 
in Table I. The case A = 0.1 and q = 0.2 is presented 
in Table II, where corresponding values of sticky tem¬ 
perature, adimensional temperature of the SW system 
and polymer packing fraction are shown. In Fig. 1 we 
present the SW potential and the effective AO potential 
that yields corresponding states for several temperatures 
and polymer packing fraction (A = 0.1 and q = 0.2). 
The inset shows the simple relation between the temper¬ 
ature of the SW particles and the corresponding packing 
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Figure 1. Corresponding states between SW potential 
(dashed lines) and AO effective potential (full lines) for A = 
0.1 and q = 0.2. From bottom to top, the inverse of the pack¬ 
ing fraction (temperature) for AO (SW) curves are 77“ 1 = 1.5, 
3 and 5.5 (T* = 0.25, 0.53 and 0.985). The inset presents the 
nearly linear relation between T* and r/ p used to map the SW 
system to the AO model. 


fraction of polymers for the AO colloid-polymer mixture 
expressed in Eq. (14). In Sec. IV, the overall approach 
will be used to analyze the spherically confined system 
of few colloids in a colloid-polymer mixture. 


III. SIMULATION METHOD 

The statistical mechanical equivalence between dif¬ 
ferent ensembles does not apply to few-body systems. 
Therefore, simulation and statistical mechanical ap¬ 
proaches should correspond to the same physical con¬ 
straints, to ensure comparable results. In this work, we 
focus on a system at constant temperature, fixed volume 
and number of particles. It corresponds to a canonical 
ensemble, and thus one assumes a Maxwell-Boltzmann 
velocity distribution. Accordingly, the simulations have 
to include a thermostating mechanism that ensures con¬ 
stant temperature and Maxwell-Boltzmann velocity dis¬ 
tribution. 

Molecular dynamics simulations of few SW confined 
in a spherical cavity were performed with a standard 
event-driven algorithm (EDMD).[53] Constant temper¬ 
ature was achieved by using a thermal-wall thermostat, 
which changes the velocity of the particle colliding with 
the wall by means of a velocity distribution compati¬ 
ble with canonical ensemble for a given temperature T. 
Thermal walls show certain features different from those 
thermostats that act over the entire system volume that 
were discussed in detail in a previous work. [12] 

For the particle interactions, we take into account dif¬ 
ferent types of events. Namely, particle-particle colli¬ 
sion and particle-wall collision. Among particle-particle 
collisions there is a further division between core and 
field events, each one related with a discontinuous step 
in <j> sw(»'). The usual EDMD algorithm was used, in 
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which the particle moves with rectilinear and constant- 
velocity dynamics, between particle collisions.[20, 53] In 
order to discriminate particle-particle collisions we used 
the logical structure of Alder and Wainwright.[54] 

The time to collision of particle i with particle j is 
calculated as: 


^ij,k — 


h 'j ± [ b ij - V ij ( r ij - a l)] 

"J, 


1/2 


(15) 


where r,y = r, — Yj and = v; — v ; are the relative 
positions and velocities of the particle pair, respectively. 
Two collision distances ox- are considered o\ = a and 
(72 = cr (1 + A). The parameter bij = r, y - • v,y must be 
negative if the particles are approaching each other and 
positive otherwise, and we consider only the positive val¬ 
ues of tij t k • The Eq. (15) is obtained by imposing the 
condition |r.y + 'VijUj^l = Ofc at collision time The 
± appears because both results are possible, given the 
right conditions: the ” sign applies for particles ap¬ 
proaching from outside regions (r.y > 02 ), and a “+” 
sign, on the other hand, is for particles already inside 
the attractive well region (r^ < ( 72 ) and getting outside 
of it, by receding. 

An ordered list of events, with increasing collision 
times tij t k is generated. Between collisions, the parti¬ 
cles move with r i = Vjt. Once a collision occurs, the new 
velocities of the pair of particles involved in the collision 
are obtained as: 


v“ ew = v° ld + Sv , 
v" ew = v° ld - <5v , 

by momentum and energy conservation dv is easily cal¬ 
culated. For core collisions 

8v = . (16) 

<7 Z 

Field cases present different possibilities. In a first 
event type, a particle enters the field reducing pair po¬ 
tential energy and consequently increasing the kinetic en¬ 
ergy: 


dv = 


2a 2 (l + \y 


^t+b 2 - 
1 11 
m J 


(17) 


Secondly, a particle is leaving the field and the pair has 
enough energy to break the bond. Then the potential 
energy increases requiring a subtraction from the kinetic 
energy: 


<5v 


2a 2 (1 + A) 2 



+ bij 


(18) 


Finally, the particles separate from each other to leave 
the field but there is not enough energy to surpass the 
well depth. Then a bounce occurs as if it were a hard 
collision 


Sv 


r ij bjj 

a 2 (1 + A) 2 ' 


(19) 


It must be considered in addition, the time at which 
each particle collides with the wall t ™. This time is cal¬ 
culated by the condition 

|U + Vitf | = R e g , 

The nearest next event is chosen as the minimum of the 
next particle and wall events: min(min(tjj ) fc), min(t^)). 
If the particle-wall collision is the next event, the system 
is evolved until the particle reaches the wall. At this 
point, the thermal-wall thermostat acts on the particle by 
imposing it a new velocity, which is chosen stochastically 
from the probability distributions: 

Pn(v n ) = m/3 |u n | exp (20) 

pt Ot) = \f^ ex p (. 

here n and t stands for the normal and tangential com¬ 
ponents of the velocities, that lie in directions — f and 
v old —(v old -f)f, respectively. The thermal walls described 
by Eq. (20) fix the temperature of the system. They 
were tested for HS confined both by planar walls and in a 
spherical pore, and produce a velocity distribution com¬ 
patible with that of Maxwell-Boltzmann[12, 55]. This 
thermal wall functionality has been extensively tested in 
a previous work with focus in confined HS particles [12] 
and we get the same results for simulations with SW par¬ 
ticles. 

Setting an N value between 2 and 6, we swep the (T, p) 
surface. For each chosen point of that surface we per¬ 
formed 10 simulations of 3 x 10® collision events, with 
a further average of the results. Temperature range was 
taken to cover two clearly distinct behavior regions. At 
low T particles rack up and form a cluster that acts as a 
rigid body, at high values particles dissociate resembling 
a HS system. Densities were picked from the very low 
“bulk like” to high values, in the vicinity of close-packing 
condition. 

From simulations we get to measure several quantities 
attained from time averages over all systems configura¬ 
tions. The studied structure and position functions are: 
the one body density function p{ r) and the averaged pair 
distribution function </(r).[56] The upper bar is just to 
distinguish our measured function from the better known 
radial distribution function g(r) that is commonly used 
to study homogeneous and isotropic systems. The cal¬ 
culated profiles for p{r) and g{r) are mean values over 
a discrete domain, obtained from a binning of spherical 
shells during the elapsed simulation time. The bin length 
is established by dividing the maximum possible distance 
value (R e g for p(r) and 2R e g for g(r)) by the number of 
desired bins, usually 1200. 


IV. RESULTS 

In this section we present the properties for the con¬ 
fined system of few short-range SW particles (A = 0.1) 



6 



Figure 2. Pair distribution function curves for N = 2 and 
p = 0.5 for different temperatures. A clear agreement be¬ 
tween theoretical and simulation results is observed. The in¬ 
set shows same curves on the range r £ [1.15, 1.97]. 

obtained from the molecular dynamics simulation. They 
are separated in structural properties, thermodynamic 
observables, and phase diagram. Structural behavior is 
analyzed in terms of the spatial correlation between par¬ 
ticles from g(r) and their spatial distribution in the cav¬ 
ity from p(r). The measured thermodynamic quantities 
describe the system as a whole, focusing on magnitudes 
that are usually utilized to characterize both bulk and 
inhomogeneous fluids composed by many-bodies. The 
statistical mechanical background was developed in Sec. 
II and in Ref. [56]. Phase diagrams should not be under¬ 
stood on the context of bulk phases. Even, they condense 
the observed system behavior in the temperature/density 
plane. For simplicity all the magnitudes are presented us¬ 
ing natural units, i.e. the unit of length is cr, the unit of 
temperature is k/e and the unit of energy is e. 

Additionally, we discuss the extent to which we expect 
an accurate mapping of different structural and thermo¬ 
dynamic properties between the AO and the SW systems. 
We analyze a few properties of the equivalent AO sys¬ 
tem based on the mapping between T and the packinf 
fraction ry p . The presented overall application of the ex¬ 
tended law of corresponding states between AO and SW 
confined system is tested at the level of phase diagram. 

A. Structural description 

We present firstly the results for the N = 2 system, 
from which exact theoretical results are available. [56] 
This allows to make a direct comparison between exact 
predictions and simulation results. By getting a perfect 
matching, we ensure that we have solid framework to ana¬ 
lyze the few-particle systems of higher N. It is also useful 
for validating the simulation program and thermostating 
procedure. The pair distribution function is shown in 
Fig. 2 for different temperatures. A clean superposition 


between theory and simulation is easily appreciable. For 
N = 2 the analytic form of g(r) is 

g(r) = Ce 12 {r) [(2 R eS - r) 2 (r + 47? e ff)] , (21) 

with ei 2 (r) = and C = 56] First, the null 

values in the range 0 < r < 1 is a trivial effect of the 
hard core repulsion. Then in the range 1 < r < 1.1 there 
is always a main peak. This peak is an expected result, 
given the shape and length of the potential: inside the 
well attractive zone, particles are more likely to be closer. 
Another common element of these curves is the presence 
of “tails”, i.e. the smooth monotonically decreasing seg¬ 
ments that are seen for ranges of r beyond the main peak. 
The tail ends at 2 R e g = 1.97 (cavity diameter) which is 
the maximum possible pair distance. Since there is not 
particle interaction for pair distance over 1.1, the tail is 
proportional to the probability of finding a pair of hard 
spheres in spherical confinement at a given distance. The 
relation between the main peak and the tail sizes is driven 
by the temperature: a steep jump in ei 2 will happen at 
low T, and ei 2 approaches to unity at high T. From a 
phenomenological standpoint, at low T the particles lack 
the energy to escape from the well, thus forming a per¬ 
manent short ranged bond. At high T, the kinetic energy 
is far greater than the well depth, meaning that the sys¬ 
tem resembles one of colliding hard cores (HS limit). For 
higher N, the competition between the structure (peaks) 
and the tail (no interaction) is one of the most visible 
effects when increasing T. The relation between g sw (r) 
and gpjoir) was not studied before, at the best of our 
knowledge. However, based on our analysis of the case 
N = 2, we propose that the corresponding functions are 
5sw(r)/ef 2 w) (r) and gxoir)/e^°\r), with r > 1. This 
mapping produces small changes in the shape of the main 
peak of gAo(r) in comparison with g sw (r). 

As can be seen in Fig. 3, where it is shown the pair 
distribution function for 3 to 6 particles inside the cavity, 
adding particles to such a small system will necessarily 
cause qualitative changes beyond the two-body analysis. 
Nonetheless, certain core aspects remain the same: those 
that are linked to the potential shape (main peak and HS 
limit). Since the integral over the complete space of g(r) 
is N(N — l)/2,[56] it is expected that lower N (at fixed 
p ) show overall lower curves. Also, R e s = (SN/Anp) 1 ^ 3 
means that at low N values and fixed number density, 
cavity size is highly susceptible to add or subtract a sin¬ 
gle particle. The shape of g(r) for different temperatures 
can be used to provide qualitative aspects of the morphol¬ 
ogy of the clusters. For very low temperatures (T < 0.1) 
the particles form a rigid cluster minimizing the overall 
system energy. Every close neighbor (1 < r < 1.1) adds 
up — e to the potential energy. The bonds are mainly per¬ 
manent, meaning stable pair distances leading to clearly 
discernible peaks. 

For these small systems, it is possible to interpret the 
low T curves g(r ) just by considering the clusters shapes. 
Three particles form a triangle and four a regular tetrahe¬ 
dron, both of which have in common that all the particles 















7 



Figure 3. g(r ) curves for TV = 2 to 6 at fixed number density 
p = 0.5. Chosen temperatures are (top to bottom in order) 
0.1, 0.25 and 0.5. The insets present the second neighbour 
peak in more detail. 


are first neighbors between themselves. This means that 
the whole probability of finding a pair at certain distance 
will localize in the first neighbor range of the pair distri¬ 
bution function (main peak). This isolated peak is shown 
for iV = 2, 3 and 4 in Fig. 3, top panel. 

The shape of g{r) for TV = 5 can be understood start¬ 
ing from the TV = 4 regular tetrahedron and then adding 
an extra particle on one of its faces. The result is an hex¬ 
ahedron composed of two regular tetrahedrons in contact 
by one face. The resulting structure has three particles, 
each one with four first neighbors and two particles with 
three first neighbors and one second neighbor. A cluster 
geometry with second neighbors gives rise to a non van¬ 
ishing probability of finding a pair distance larger than 
the main peak region. The stable structure with the sec¬ 
ond neighbor distance in a constrained region results in 
a second peak. 

For the case N = 6 there are two observed cluster ge¬ 
ometries. One of them is the regular octahedron. How¬ 
ever, the most frequent geometry observed in the simu¬ 
lations is an irregular octahedron. This irregular poly¬ 
hedron has a typical path of formation starting from a 
hexahedral cluster of five particles to which the remain¬ 
ing particle adds over one of its faces. For extremely low 


temperatures (T <C 0 . 1 ), once the particles form their 
bonds, they will stay bonded permanently. For a softer 
cluster (T < 0.1), single bonds have a slight chance of 
breaking and a rearranging of the cluster structure can 
take place. The secondary peak of g[r) is sensitive to 
these different structures as shown in the inset of the top 
panel in Fig. 3. There are two secondary peaks both 
related with second neighbor distance: the larger one 
represents the second neighbors in the irregular octahe¬ 
dron while the smaller one is characteristic of the regular 
body. The significant height difference shows that the 
irregular cluster is more frequent over time. We point 
out also that both geometries have the same number of 
first neighbors being thus isoenergetic. Therefore, from a 
statistical mechanics point of view, the only factor that 
can lead to the prevalence of one geometry over the other 
is strictly coming from entropic contributions. The irreg¬ 
ular cluster presents lower symmetry and higher entropy. 

Increasing the temperature leads to the break up of 
multiple bonds, which results in flexible or plastic clus¬ 
ters. For intermediate T ranges (0.1 < T < 0.25), the 
particles remain constantly linked, but now they are not 
tightly bound. Certain bonds are likely to break, en¬ 
abling the particles to displace inside the cluster. The 
inset of Fig. 3 for T = 0.25, shows that the isolated 
peaks are surrounded by non-vanishing values. For ex¬ 
ample, the triangle for TV = 3 breaks one of its pair bonds 
in such a way that it opens up and stretches like a chain. 
For higher TV, is essentially the same. Additional transla¬ 
tional freedom leads to possible deformations into wider 
shapes than the original rigid body. The weaker the 
structure is, the lower the main and second peak become, 
favoring pair distance probability on the surrounding ar¬ 
eas. 

For higher temperatures, (T > 0.3) soft cluster starts 
to dissociate and single particles are free to have any 
distance from the cluster, inside the cavity limits. Any 
kind of stable structure fades out and only instantaneous 
pairs endure. As can be seen in the bottom panel of Fig. 
3, the tail engulfs any close range structure and the main 
peak reduces around half of its height, compared to the 
T = 0.25 panel. Higher temperatures do not add any 
qualitative variation: the main peak will decrease until 
it becomes part of the tail, in the HS limit. 

Density profiles, shown in Fig. 4, present very explic¬ 
itly the inhomegeneity of the system. Unlike commonly 
studied bulk systems, there are significant local spatial 
variations of the one body density function. Note that 
R e s varies with TV, in Fig. 4. It takes values from R e g = 1 
for two particles to R c g = 1.9 for six particles. For the 
case TV = 2 the density vanishes in the center because 
if one particle in placed at r ~ 0 , the available volume 
for the other one becomes very small. For density val¬ 
ues that define effective radii higher than the particle 
diameter, i.e. TV = 3 to 6 , all the profiles have simi¬ 
lar form, independently of particle number. The profiles 
have two distinct regions. An approximately constant 
density on the center of the cavity, that is cut at the 































Figure 4. Density profiles for N = 2 to 6 at fixed number 
density p = 0.5. The chosen temperatures are (top to bottom 
in order) 0.1, 0.25 and 0.5. 

vicinity of the wall (R e g — r < 1) and the "interfacial" 
region closer to the wall. The extension of the plateau 
depends on the relation of cavity size and particle diam¬ 
eter, as observed from the different sizes at equal number 
density in Fig. 4. For low temperatures the system may 
be treated as a single nearly-rigid body, free to translate 
and rotate in the central region. When the cluster gets 
closer to the wall, some possible cluster orientations are 
restricted, leading to a reduction in rotational entropy. 
Consequently, it is more likely to find the cluster in the 
central region. Higher temperatures soften the particles’ 
bonds progressively, causing a reduction of the dispar¬ 
ity between the plateau and the region close to the wall. 
Once dissociation becomes dominant for high tempera¬ 
tures, depletion arises and the wall starts to have a more 
intense effective attraction. Further increase in the tem¬ 
perature makes the wall attraction higher, while the par¬ 
ticle correlation disappears. At the HS limit the highest 
probability of finding a particle is at the wall. This was 
also observed in a system of pure HS particles in spheri¬ 
cal confinement.[12] For the AO system we expect similar 
density profiles to those of SW for equivalent tempera¬ 
ture and packing fraction. The temperatures T = 0.1, 
0.25 and 0.5 correspond to (top to bottom panels in Fig. 
4) polymers inverse packing fraction rj= 0.694, 1.50 
and 2.84, respectively. As we will see, the case T = 0.1, 
may be a too small temperature to use the extended law 
of corresponding states. 

The structure dependency on p is far more subtle than 


Figure 5. g(r) curves for N = 2 to 6 at fixed temperature 
T = 0.25. Number densities are 0.1 (Top) and 0.8 (Bottom). 
The insets present the second neighbour peak in more detail. 


on T. In Fig. 5 curves of g(r) at T = 0.25 and for two dif¬ 
ferent densities are shown. One observes that the shape 
of the peaks remain practically unaltered. Global values 
of the curve raise for higher number densities, because 
the normalization is the same in a smaller cavity size. 
We have shown the case T = 0.25, as an example but we 
observed the same general picture for other values of T 
which are not presented here. 

From here on, we select the case N = 5 to give more 
detailed analysis. This case is the only one that presents 
a second peak but does not have multiple stable geome¬ 
tries, as the case N = 6. At low temperatures it is ob¬ 
served a clearly defined structure, while retaining little 
longer range order. The N = 5 results can be extrapo¬ 
lated to the other few-particle systems. Density profiles 
are shown in Fig. 6 for a wide range of temperatures 
and number densities. As already pointed out, at low 
and intermediate densities there exist a plateau and an 
interfacial region close to the wall. Increasing the temper¬ 
ature leads to an overall probability density favoring po¬ 
sition closer to the wall, as a result of relatively stronger 
depletion attraction. [12, 36] The higher the density, the 
smaller the available free volume for the rigid cluster, 
thus the plateau becomes smaller and steeper. At cer¬ 
tain number density there is no more room to locate a 
particle in the center of the cavity. The available space 
is so small that a particle at the center would push the 
remaining ones out of bounds. This is shown in the lower 
panel of Fig. 6, which exhibits an excluded volume region 
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Figure 6. Density profiles for TV = 5 at different tempera¬ 
tures. Number densities are 0.1 (Top), 0.8 (Center) and 1.1 
(Bottom). 


close to the center of the cavity. This gives a clear insight 
of how, at high densities, the system conformation and 
translation becomes dominated by the cavity shape. At 
high densities, rigid clusters have low translational free¬ 
dom and therefore their constituent particles stand at 
approximately fixed distances from the center. Cluster 
structure gets expressed on the density profiles that show 
local maxima and minima. Increasing the temperature 
softens the cluster, causing the local structure features to 
disappear, leaning towards monotonous curves. Finally, 
at dissociation temperatures, the depletion dominance 
gets clear and the density at the wall is the highest in 
the profile. We note here a difficulty, intrinsically related 
with the spherical shape of the cavity. The local proper¬ 
ties are hard to measure in the central region because it 
is poorly sampled. Indeed, the sampling becomes poorer 
as r decreases towards the center of the pore. This ef¬ 
fect is produced by the rapid reduction of the sampled 
volume that produce large fluctuations in p(r) and other 
quantities. These fluctuations can also be observed in 
Fig. 4. 

In Figure 7 a similar systematic approach is applied 
for the g(r) curves under density variation. Each studied 
density p = 0.1, 0.8 and 1.1 corresponds to a cavity di¬ 
ameter of 4.57, 2.28 and 2.05, respectively. By increasing 
the temperature, pair bonds are more likely to break and 
produce dissociations. The odds of a separated pair to 
become together again is smaller at larger free space. As 
a consequence, main peaks fall more abruptly at dissocia- 



Figure 7. g(r) for N = 5 at different temperatures and rep¬ 
resentative number densities. Left panel: detail of the main 
peak, right panel: detail of the second peak for r € [1.15; 3]. 


tion temperatures and lower densities. At high densities, 
there is not enough room for the particles to stay away 
from each other, which forces an increase in the prob¬ 
ability of finding pair separations inside the well range, 
even if particles have a very high kinetic energy, as com¬ 
pared to the well interaction energy e. At p = 1.1 it is 
noticeable how cavity diameter is close to second neigh¬ 
bor distance. For N = 5 the rigid cluster geometry is 
also the most compact the system can achieve, the cav¬ 
ity is only a slightly larger than the smallest possible 
configuration. This implies that the cluster as a whole is 
practically locked in the center, allowed mostly only to 
rotate. As already observed in Figure 6, at high densities 
and low temperatures, the particles of the rotating rigid 
body maintain a stable distance from the center. For 
p > 1.1, the wall cavity squeezes the system, shortening 
the second neighbor distance. 

Despite the SW interaction is isoenergetic once inside 
the field range, the main peak of g(r) shows a negative 
slope as can be observed in Figs. 2, 3, 5 and 7. This 
implies that, within the interaction range r E [1,1 + A], 
particles tend to be at closest distance instead of near the 
external side of the well. Theoretical results for N = 2 
(Eq. 21) points out that the slope originates from the 
particle-cavity interaction. Basically the main peak is 
an offset mounted on a decreasing function. This can 
be rationalized by noting that a closer pair has more 
free space in the cavity than a stretched one, increasing 
the traslational entropy. Also particle collisions with a 
curved concave wall will, in average, tend to group them 
together. 


B. Thermodynamic quantities 

We focus on four thermodynamic properties to char¬ 
acterize the system as a whole. We analyze the pressure 
and the energy of the system, that have robust definitions 
and are measured in a straightforward way. Additionally, 
we study the surface tension and the surface adsorption, 
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Figure 8. Compressibility factor /3P w /p as a function of tem¬ 
perature, for different number densities for N = 3 (Top panel) 
and N = 6 (Bottom panel). 


which are intrinsically related with the inhomogeneous 
nature of the system. These last quantities are more 
subtle and difficult to measure by simulation. 

Pressure on the wall is expected to increase with tem¬ 
perature and number density, given that both parameters 
should increase the average number of collisions on the 
wall. We work with the compresibility factor (3P w /p in 
Fig. 8 to eliminate the linear dependence, allowing to dis¬ 
tinguish deviations from the ideal case /3P W /p = 1. Only 
the high temperature and very low density cases follow 
the ideal behavior, when the relative well attraction is too 
weak and excluded volume from the cores is negligible. 
Cluster to dissociation temperatures are mediated by a 
sudden increase in f3P w / p and further saturation. The 
jump becomes smoother for higher densities as a result 
of what has been pointed out from g(r) analysis: with 
less space for separation, qualitative differences between 
a cluster and unbounded particles are smaller. The men¬ 
tioned depletion emergence, for any density value at high 
T, explains the asymptotic increase of [3P w /p until the 
HS limit. 

At very low temperatures, the system is far away from 
the ideal gas behavior because the relative potential well 
produces a cluster. In this limit the system compress¬ 
ibility factor increases monotonically for increasing den¬ 
sity. We attribute this to the behavior of a unique 
finite-size cluster allowed to stay in an effective volume 
^cluster < V. The single cluster behaves as an ideal gas, 
and thus, its compressibility factor is /3.PI4iuster = 1 
i.e. flP/p — IV _1 Vy Wiuster- This explains the low- 



Figure 9. E — E c0 h as a function of T for different number 
densities Top panel shows N = 3 and Bottom panel N = 6 
cases. The horizontal dashed line plots |S co h| = —E co h. 


density and low-temperature values /3P w /p > 0.33 and 
PP W / P > 0.17 for N = 3 and 6 respectively. At interme¬ 
diates temperatures, when the bond-breaking probability 
is non-negligible, the case of low density raises much more 
rapidly to its saturation value than those of higher den¬ 
sities. This could be related to the strong reduction of 
recombination rate at lower densities. Once a bond is 
broken the probability of the particles to meet again is 
very small. This is not the case for intermediate to high 
densities. 

The mean energy per particle E is the addition of 
the kinetic term 3T/2 driven by the temperature and 
the mean potential energy. For low temperatures, where 
\4>n\/N 3T/2, the system is in a rigid cluster state. 

Then, E can be precisely calculated as the number of 
bonded pairs for a given geometry. We call that value 
cohesion energy E co h, which is the lowest (fundamental) 
possible energy of the system and we define it as the zero 
value in Fig. 9, presenting the energy versus tempera¬ 
ture. Average energy per particle has a similar behavior 
to the one of the pressure in Fig. 8. For low density, it 
presents an abrupt increment in going from low to higher 
temperatures (T ~ 0.2). This jump agrees with the range 
of non-rigid cluster, ending at dissociation temperatures 
(T ~ 0.3). Then it follows a weak linear variation for 
high temperatures, according to the equipartition theo¬ 
rem. It is worth noting that for systems with a unique 
stable rigid cluster geometry at low temperatures, all the 
curves must collapse to a single with slope 3/2, inde- 































Figure 10. Adsorption T as a function of temperature for 
different number densities. The cases N = 3 (Top panel) and 
N = 6 (Bottom panel) are shown. 

pendently of the density. For reference, the behavior of 
both, the pure kinetic energy and the shifted one (plus 
| f? co h I); are also shown in dotted and dot-dashed lines 
in Fig. 9. Rigid clusters constitute compact structures, 
and increasing the density does not changes the number 
of first neighbors. Higher density lines have lower values 
because the particles are closer, so SW interactions are 
forced. The changes produced with increasing tempera¬ 
ture in each curve are more pronounced for higher TV, be¬ 
cause the rigid cluster has more bounds (per particle) to 
be broken. This feature is also shown for the different val¬ 
ues that takes |.E C oh|> hr dashed line. The line also serves 
to visualize the condition E = 0, which characterize the 
equilibrium between potential and kinetic energies. This 
crossover line separates two characteristic regions: one 
where the potential energy dominates over the kinetic 
energy, proper of clusters, and another one where kinetic 
energy dominates, a feature proper of gases. 

Surface adsorption F and surface tension 7 are ba¬ 
sic properties used to characterize the inhomegeneity in¬ 
duced on the system by the presence of walls. They are 
measured in the same way as a in a previous work. [12] We 
won’t delve into details and only give here the definition 
of F, and the expression of 7 , based on Laplace equation: 
T = (p- Pc )E and 7 = Pc ~ Pm R e g. Here p c (P c ) refers to 
the average density (pressure) near the center of the cav¬ 
ity. These magnitudes are difficult to measure because 
one must fix a criterion to choose the region where aver¬ 
ages should be done. The criterion must be applied for 
all the available range of p and T. Note that at large R e g, 



Figure 11. Surface tension as a function of T for different 
number densities. The cases N = 3 (Top panel) and N = 6 
(Bottom panel) are shown. 


the density profiles attain a nearly constant value in the 
central region (plateau). At constant N, the density p c 
becomes higher for smaller cavity radius until a limiting 
value, in which particles cannot freely place themselves 
in the center. This central region progressively turns into 
an excluded zone. It is important to note that at such 
high packing values, the definition of p c is less accurate, 
since there is not a clear distinction between the center 
region and that in the vicinity of the wall. Adsorption is 
shown in Fig. 10. For low densities, adsorption is neg¬ 
ative for small T and then experiences a jump around 
dissociation temperatures to flatten at higher values of 
T. It becomes positive at high temperatures. For high 
densities, adsorption is positive and nearly independent 
of the temperature. The limiting cases of temperature 
are clearly identified. For high temperature F is positive 
and increase monotonously with density. In the case of 
low temperature, for small densities T is negative and de¬ 
creases with increasing p, up to a certain minimum value. 
Further increment of p produces a sudden rise of T, that 
becomes positive. This behavior goes in line with that 
observed in the density profiles, in Section IV A. p(r) 
presents an enhancement close to the wall for F > 0 and 
an increment in the center for T < 0 . 

Fig. 11 shows the surface tension 7 for three (Top 
panel) and six particles (Bottom panel). At vanishing 
temperature they start from 0 , having then two distinc¬ 
tive behaviors. For low to intermediate densities the 
7 curves are positive at low T. They start with posi¬ 
tive slope, reach a maximum value, to become negative 
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at high temperatures. For high densities, 7 curves are 
negative. They start with negative slope and decrease 
monotonously with temperature. 

The two different characteristics of 7 can be rational¬ 
ized by considering that when clusters are favored at 
medium to low densities, the system tends to get far from 
the cavity wall, having the cluster size as its characteris¬ 
tic size. This minimizes the intrinsic area of the system, 
going along with a negative adsorption and an increase of 
density at the center of the cavity. At high enough tem¬ 
peratures, the system behaves as HS particles, having the 
confining cavity as a characteristic size, with an effective 
entropic attraction from the wall, and a negative surface 
tension. [12] The curve that shows the global maximum of 
surface tension corresponds to higher density for N = 6 
than for N = 3. Also, the temperature of those global 
maxima of 7 is shifted towards higher values. 

The approximate linear dependence of 7 with T for 
both very low and high temperatures is explained by the 
expected hard sphere limit where /3y only depends on 
density. The behavior upon variation of density is similar 
to that observed for T. 

Some of the measured thermodynamic properties of the 
SW system can be readily mapped to the AO model by 
the established relation between T and ?y p . We expect one 
of these magnitudes to be the pressure on the wall, once 
the energy scale is compensated [as in the case of (3(j) ao 
in Eq. (9)], thus j3P w but also f3P w /p that was plotted in 
Fig. 8 could be mapped. The measured surface tension 
is similar to the pressure, it should be transformed to 
/Fy. A third magnitude is T which does not scales with T 
and depends on characteristic features of p{r). The case 
of energy is more complicated because in the AO system 
the energy is purely kinetic, and therefore the energy can 
not be mapped. 


C. Phase diagrams 

In Fig. 12 we show the change of the main charac¬ 
teristics of the phase diagram with the variation of N. 
This summarizes variation of structural properties with 
T and p. The vertical lines come from the analysis of 
the pair distribution function, by comparing the height 
of the main peak g(r = 1 + ) with that in the region next 
to the main peak g(jr = 1.1 + ), and its connection with 
the qualitative behavior of the system, obtained by di¬ 
rect visualization of the dynamics of the particles. C 
corresponds to g{r = 1.1 + )/g{r = 1 + ) = 0.005 and D 
to g(r = 1.1 + )/g(r = 1 + ) = 0.03. The relations between 
those points have been picked by noting that those values 
match cluster softening (C) and dissociation process (D). 
These lines divide temperature domains by particle con¬ 
formation: low temperatures until C mostly defined by 
particles gathering in a hard cluster. Between C and D 
the system forms a soft cluster with relative movements 
among particles. For higher temperatures, beyond D, 
frequent dissociations are observed. The horizontal lines 
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Figure 12. Superimposed phase diagram for N = 2 to 6. Ver¬ 
tical lines come from analyzing the pair distribution function. 
At the left side of C, there is a hard cluster region. Between 
C and D we define a soft or plastic cluster region and at the 
right side of C dissociation starts to happen. Horizontal lines 
show density points where density profiles are maximum (B), 
or vanish at the center (A). The rectangles show the range of 
p between A and B positions for every N. 

come from the analysis of the main features of the den¬ 
sity distribution. Line A indicates the cases in which 
the density profile at the center of the cavity p(r = 0) 
reaches its maximum, and line B when it becomes zero 
\p{r = 0) = 0]. The former case represents low trans¬ 
lational freedom for the particle that is at the center, 
while for the later there is such a high density (or a small 
cavity), that a particle can not be in the center due to 
the excluded volume. These horizontal lines delimit the 
following regions: below B there is a zone in which the 
system is moderately inhomogeneous. The region limited 
by A and B corresponds to a strong reduction of freedom 
of motion that reduces the occupation in the center. For 
densities beyond line A the center of the cavity is an 
excluded volume (high confinement). Note that a third 
horizontal line (not shown) fix the maximum density of 
the confined system where it becomes completely caged. 
This density can be calculated with a simple geometrical 
approach and varies with N. In contrast to the vertical 
lines, the position of A and B are strongly dependent on 
N. 

In Figs. 13 to 17 we present the phase diagram for the 
different number of particles studied. It is shown a set of 
characteristic values of density and temperature labeled 
with letter A, B, C and D. They separate regions where 
the systems present different qualitative behavior. These 
regions would represent different phases of the system in 
a macroscopic context. 

Additionally, the scalars chosen to distinguish different 
qualitative behavior of the system are the (T, p) points on 
which: kinetic energy is equal to potential energy (black 
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Figure 13. Phase diagram for N = 2 with characteristic 
curves. A, B, C and D thick lines are the same from Fig. 12. 
Additionally, curves indicating E = 0 (black circles), f3P w = p 
(open red squares); T = 0 (green diamonds), and 7 = 0 (blue 
triangles) are shown. Red crosses present the values obtained 
by interpolating the curves of j3P w / p for specific density val¬ 
ues. These lie between points for which simulation results are 
available thus allowing to complete the curve at the saturation 
area. The two red filled squares show exact analytic values for 
f3Pw = p and T = 0 at p —► 0. Green bigger squares represent 
the AO model results for f3P w = p. In the top x-axis the scale 
of packing fraction of polymers in the AO model is presented. 


circles), pressure on the wall follows ideal gas equation of 
state (red open squares and crosses), adsorption and sur¬ 
face tension change their signs (green diamonds and blue 
triangles, respectively). We point out that the (T,p) grid 
step is 0.1, except for the cases of low temperatures (e.g. 
T = 0.01) and higher densities where the steps become 
more spaced. In particular, low temperature f3P w = p 
points (red crosses) have been attained by interpolating 
/ 3P W curves for specific density values, whose variations 
are too small for our scale. 

The E = 0 curve (open circles) divides the phase space. 
Towards the left, the (modulus of) potential energy is 
higher than the kinetic energy and the opposite case is 
true, towards the right. In the region with E < 0 (at 
the left of E = 0) the potential energy dominates. This 
feature is characteristic of systems that spontaneously 
collapse in cluster aggregates or condensed phases, but 
also of cold enough compressed systems, where the avail¬ 
able space is reduced. On the opposite, when E > 0 (at 
the right of E = 0 curve) the kinetic energy dominates 
and this is typical of systems of free particles, as the case 
of diluted gases. Note that this last region includes a 
high enough temperature and high density range, where 
particles stay together, as a consequence of the strong 
confinement. The E = 0 curve is specially sensitive to 
changes in N, since it is affected by the average number of 
bonds per particle. At low IV, adding particles increases 


Figure 14. Phase diagram for N = 3 with characteristic 
curves. Symbols, lines and top x-axis follow the same no¬ 
tation from those of Fig.13. 
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Figure 15. Phase diagram for N = 4 with characteristic 
curves. Symbols, lines and top x-axis follow the same nota¬ 
tion from those of Fig. 13. The inset shows details of /3P W = p 
at low temperatures. 


the maximum number of bonds per particle, given that 
the system is too far away from bulk case. As a result of 
this, the curve shifts towards higher temperatures upon 
increasing N. This is specially observed comparing Fig¬ 
ures 13 and 14. The curve flattens for higher densities 
since smaller cavity forces bond formation, leading to a 
fixed potential energy value of a tightly packed cluster. 
This curve has not a counterpart for the colloid-polymer 
mixture analyzed with the AO model. Even when there 
is an effective potential between colloids, the system is 
athermal and the origin of the potential is purely en- 
tropic. Both, colloids and polymers have only kinetic 
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Figure 16. Phase diagram for N = 5 with characteristic 
curves. Symbols, lines and top x-axis follow are the same 
notation of those of Fig. 13. 
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Figure 17. Phase diagram for N = 6 with characteristic 
curves. Symbols, lines and top x-axis follow the notation of 
Fig. 13. 


energy, thus E yZ 0 for non zero temperatures. 

The curve (3P W = p (open squares) indicates the points 
of the (T, p) where the EOS for the system pressure be¬ 
haves like the ideal gas pressure T\,\. At the left side and 
below of the 0P W = p curve, P w < Pid, the system is un¬ 
dercompressed with respect to the ideal gas (at the same 
temperature and number density). This region encloses 
the origin. Outside this region P w > P id the system is 
overcompressed in comparison with the ideal gas. For 
low densities and high temperatures this is expected as 
a consequence that depletion favors wall contact. Higher 
densities force wall contact, independently of tempera¬ 
ture. The theoretical prediction based on the limit of 


large V (low density) is shown with a red square on T 
axis. Notably, the obtained limiting temperature fits the 
curve and is independent of N. From N = 3 to 6 all these 
curves follow a similar behavior: from zero density up to 
p w 1, where T ~ 0.3. The curves at lower tempera¬ 
tures show a strong dependence with N. It is interesting 
to note that the critical temperature of the vapor-liquid 
metastable transition in the studied short-range SW sys¬ 
tem is T = 0.47, which, for the largest values of N, nearly 
coincides with the intersection of zero-energy and ideal 
gas pressure curves. [57] 

The 6P W = p curve was selected to verify the validity 
of the corresponding-states mapping between SW and 
AO in confined systems. Using Metropolis-Rosembluth 
Monte Carlo calculations[58, 59] we have evaluated the 
density distribution of the AO system (for few values of 
p), and used the contact theorem Eq. (11) to evaluate 
j3P w . Thus, we seek for the value of r/ p that produce 
/3P W = p. Extended corresponding state law [see Eq. 
(14)] was used to evaluate the temperature of the corre¬ 
sponding SW equivalent system. Calculated values are 
drawn in green squares at Figs. 13 to 17. The obtained 
values of 77“ 1 are given at the top horizontal x axis, while 
the equivalent temperature of the SW system can be read 
at the bottom axis. We found a general coincidence be¬ 
tween Monte Carlo results for AO and simulation results 
for SW, along all the analyzed range. At lower temper¬ 
atures, near clusterization, the mapping between both 
systems becomes poorer. This is expected, because the 
use of the extended law of corresponding states is not 
documented for freezing temperatures. 

Zero surface adsorption and surface tension curves in¬ 
dicate where there is no excess in surface of particle con¬ 
centration or free energy, respectively. As it was men¬ 
tioned in Sec. IV B, both are difficult to measure. T 
and 7 reveal a strong reduction of accuracy at very low 
density where the system is quasi-homogeneous and both 
quantities become very small. This makes specially hard 
to measure the value of T at which r = 0 and 7 = 0 for 
p —> 0. In addition, given that T and 7 become ill defined 
at high densities we do not evaluate zero surface excess 
in this case. In summary, the results presented here give 
a general idea of the position and shape of both curves in 
the (T, p) plane. Adsorption illustrates the relation be¬ 
tween the density of particles at the center of the cavity 
and those closer to the wall. Towards the left side, where 
T < 0, the particles are more likely to be at the center 
of the cavity. For T > 0, at the right side, particles favor 
positions close to the wall. 


V. CONCLUSIONS 

In this work we studied thoroughly the properties of 
few colloidal particles confined in a spherical cavity. We 
adopt the short-range square well model and provide a 
deep characterization of the structural and thermal prop¬ 
erties of systems of 2, 3, 4, 5 and 6 particles in a spherical 
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pore for the complete relevant ranges of density and tem¬ 
perature. Additionally, we compare the simulations with 
exact results for the case N = 2. Applying an extension 
of the corresponding states law to confined systems, we 
establish a mapping between the square well system and 
the AO model for effective interactions between colloids 
in a polymer colloid mixture. We also developed the sta¬ 
tistical mechanical approach to systems of few particles, 
SW and AO, in confinement and map the pressure on 
the wall at a given temperature in the SW system to 
the equivalent packing fraction of polymers in the AO 
system. 

The structure of the system ranging from low density 
to almost caging of the particles in the cavity was char¬ 
acterized through the pair correlation function and den¬ 
sity profiles for the entire relevant range of temperatures. 
Thermal bulk properties such as energy and pressure on 
the wall were calculated and characterized. Different ef¬ 
fects of confinement were also studied, identifying their 
energetic or entropic origin and focusing on the inhomo¬ 
geneities present in the system. Surface properties were 
analyzed with quantities reminiscent of surface tension 
and adsorption in macroscopic counterparts of the square 
well system. 

We characterized the morphology of these systems, 
defining different regions of similar behavior and crite¬ 
ria to provide phase diagrams in the (T, p) plane, for the 
different number of particles. In this phase diagram, we 
identified temperature regions where the system behaves 
as a rigid cluster, as a plastic cluster, and a region where 
the system dissociates, up to the limit of hard-sphere- 
like behavior at very high temperatures. In the density 
domain, we recognized regions with different degrees of 
inhomogeneity which can be classified in the following 
categories: low-to-moderate, moderate to excluded vol¬ 
ume, and excluded-volume to caging regions. We defined 
several characteristic curves in the phase diagram, such 
as that of zero energy, ideal gas pressure, zero adsorption 
and zero surface tension. These lines delimit meaning¬ 
ful references in the phase diagram, that were used as a 
complement for the analysis. 
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are the spherically symmetric pair potentials. Note that 
the region C where the center of colloids are confined 
(with volume V ) is different to the region V where the 
center of polymers lies (with volume V p ). In fact, the 
boundary of V must be placed in a region where the 
polymers reach their bulk properties. For the AO system 
with q < 0.1547 confined in a spherical pore, the smallest 
region V is an sphere with radius R$ + a p /2. 

Polymers behave as ideal gas particles. If we fix the po¬ 
sition of the colloids, they exert a fixed external potential 
to the polymers and thus 


rNp 


e 04>cp drp p = 
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N v 


= z, 


N„ 


(A2) 


where Z® is the Cl of one polymer in V p at fixed colloids. 
Furthermore, 




, (A3) 


= [ e {-d4>c.+z P z^) dr N , ( A4 ) 

Jv N 


and thus we can simply analyze the case of one polymer. 
We introduce the Mayer function for the colloid/polymer 
Boltzmann statistical weight e* = exp(—= 1 + /* 
in Z ( g, to obtain 


Z <s = (l + X! fifi + --)dr p , (A5) 

^ Vp i <ij> 


where higher order terms are products of three or more 
functions / concerning the position of three or more col¬ 
loids. In this integrand /,; is minus one for r p such that 
the (center-to-center) ith-colloid to polymer distance ful¬ 
fills Xi < a + (J p and otherwise is zero, f jj'j is one if r p 
fulfills both Xi < a + <jp and Xj < a + a p and is zero oth¬ 
erwise, and so on. Once integrated, Z t g> takes the form 


Financial support through grants PICT-2011-1887, 
PICT-2011-1217, PIP 112-200801-00403, INN-CNEA 
2011, PICT-E 2014, is gratefully acknowledged. 


Appendix A: Effective potential and AO partition 
function 


^ ^ y (T 3 (1 + qf N + J2 V 0 (r«) + ... (A6) 

<ij> 

Here V 0 (rij ) is the overlap volume between two spheres 
with radius (cr + cr p )/2 and extra terms include the over¬ 
lap of at least three spheres. Turning to the integrand of 
Eq. (A4), we utilize the identities z p = p p and z p V p = N p 
to obtain 


We show here how to transform Eq. (7) in Eq. (8). To 

z Np 

this end, we analyze the term ^2 N ~^-\Zn,n p in Eq. (7), 
where the (canonical ensemble) colloid-polymer mixture 


exp (N p - N^) exp -/3 ^ <fes (nj) + P P ^1 Vo(nj) , 
<ij> <ij> 

(A7) 



16 


with N* = PpV e xcN and i> exc = ^cr 3 (1 + q) 3 ■ In addi¬ 
tion, we neglected higher order terms in Eq. (A7). These 
terms are null if q < 0.1547. For q > 0.1547, including the 
case q = 0.2 analyzed in the present work, one expects 
that three-body contribution will be negligible in com¬ 
parison with two-body terms. Naturally, 'B.p.h = exp N p 
and /34>Ao{r) = /3(pHs{r) - p p V 0 (r), where /3</>Ao(r) is the 
same expression given in Eq. (9). Therefore, Eq. (A4) is 

Z^ 0 ' 1 . This demonstrates the equivalence between Eq. 
(7) and Eq. (8). 

The described procedure can be generalized in several 


ways. It is not restricted to the spherical pore, and thus, 
it applies to other pore geometries like cylinders, cuboids, 
slits, single walls, etc. Further extensions include the 
case of non-free polymers where both, colloids and poly¬ 
mers, are confined, the AO model with q > 0.1547, and 
also other non-AO systems with more general interac¬ 
tion potentials. It can also be readily applied to sys¬ 
tems of AO particles in spaces with dimensions other 
than three (discs and hyper-spheres), with V (> (r) taken 
from Ref. [60]. 
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